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Abstract 

We perform a bifurcation analysis of the mathematical model of Jones &: 
Kompala [K.D. Jones and D.S. Kompala, Cybernetic model of the growth 
dynamics of Saccharomyces cerevisiae in batch and continuous cultures, J. 
Biotech., 71:105-131, 1999]. Stable oscillations arise via Andronov-Hopf bi- 
furcations and exist for intermediate values of the dilution rate as has been 
noted from experiments previously. A variety of discontinuity induced bifur- 
cations arise from a lack of global differentiability. We identify and classify 
discontinuous bifurcations including several codimension-two scenarios. Bi- 
furcation diagrams are explained by a general unfolding of these singularities. 



1 Introduction 

Yeasts are single-celled fungi of which more than one thousand different species have 
been identified. The most commonly used yeast is Saccharomyces cerevisiae which 
has been utilized for the production of bread, wine and beer for thousands of years. 
Biologists in a wide variety of fields use S. cerevisiae as a model organism. 

A common experimental method for observing biochemical processes involved in 
yeast growth is that of continuous cultivation in a chemostat [1]. The cell growth 
takes place in a vessel that is continuously stirred. A nutrient containing fluid 
is pumped into the vessel and cell culture flows out of the vessel at the same rate, 
ensuring that the volume of culture in the reaction vessel remains constant. The rate 
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of flow in and out divided by culture volume is called the dilution rate. Quantities 
such as concentrations of chemicals can be measured in a variety of ways, see [2], [3] 
for methods used in S. cerevisiae experiments. 

As a continuous culture experiment is carried out, it is common for the system 
to reach a steady state. At the steady state, the rate of cell division in the culture 
is equal to the dilution rate. However experimentalists in the late 1960's, [4], ob- 
served that instead of settling to a steady state, continuous culture experiments of 
S. cerevisiae could in some cases produce stable oscillations. Von Meyenburg, [5], 
discovered in subsequent experiments that these oscillations only occur in an inter- 
mediate range of values of the dilution rate (between about 0.08h~^ and 0.22/i~^). 
Much work has since been done to understand the cause of such oscillations, see for 
instance PIUE]. 

S. cerevisiae has three metabolic pathways for glucose: fermentation, ethanol 
oxidation and glucose oxidation. The model of Jones & Kompala [9] hypothesizes 
that the competing metabolic pathways of the growing yeast cells create feedback 
responses that produce stable oscillations. It assumes that micro-organisms will 
utilize the available substrates in a manner that maximizes their growth rate at 
all times. To enforce this optimization a "maximum function" is introduced in 
the model equations; as a result, the model is an example of a piecewise-smooth, 
continuous dynamical system. 

Piecewise-smooth systems are characterized by the presence of co dimension-one 
phase-space boundaries, called switching manifolds, on which smoothness is lost. 
Such systems have been utilized in diverse fields to model non-smooth behavior, for 
example vibro- impacting systems and systems with friction p!0l[Tn[T2] . switching in 
electrical circuits [l3l|lll[T5], economics [Ml ITT] and biology and physiology, p^ [T9]. 

The interaction of invariant sets with switching manifolds often produces bifur- 
cations that are forbidden in smooth systems. For instance, though period-doubling 
cascades are a common mechanism for the transition to chaos in smooth systems, 
in piecewise-smooth systems periodic orbits may undergo direct transitions to chaos 
pm [21] . These so-called discontinuity induced bifurcations can be nonsmooth ana- 
logues of familiar smooth bifurcations or can be novel bifurcations and unique to 
piecewise-smooth systems. A bifurcation in the latter category that is simple in 
appearance (for example the transition from a stable period-1 solution to a stable 
period-3 solution in a piecewise-smooth map) often corresponds to a combination 
of or countable sequence of smooth bifurcations. In this situation, arguably, the 
piecewise-smooth system describes the dynamics more succinctly than any smooth 
system is able to. Alternatively bifurcations in piecewise-smooth systems may be 
extremely complicated, see for instance [201 [121 HI] and references within. 

A piecewise-smooth, continuous system is one that is everywhere continuous 
but nondifferentiable on switching manifolds. In such a system, the collision of a 
mathematical equilibrium (i.e. steady state, abbreviated to equilibrium throughout 
this paper) with a switching manifold may give rise to a discontinuous bifurcation. As 
the equilibrium crosses the switching manifold, its associated eigenvalues generically 
change discontinuously. This may produce a stability change and bifurcation. In 
two-dimensional systems, all co dimension-one discontinuous bifurcations have been 
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classified [22j , but in higher dimensions there are more allowable geometries and no 
general classification is known. See for instance [231 El] for recent investigations into 
three-dimensional systems. 

In this paper we present an analysis of discontinuity induced bifurcations in the 
eight- dimensional S. cerevisiae model of Jones & Kompala [9j. The model equa- 
tions are stated in ^ In ^we illustrate a two-parameter bifurcation set indicating 
parameter values at which stable oscillations occur. The bifurcation set also shows 
curves corresponding to co dimension-one discontinuous bifurcations. These bifur- 
cations are analogous to saddle-node and Andronov-Hopf bifurcations in smooth 
systems. Bifurcations relating to stable oscillations are described in ^ We observe 
period-adding sequences over small regions in parameter space. In ^ we provide 
rigorous unfoldings of codimension-two scenarios seen in the bifurcation set from a 
general viewpoint. Finally [|6] presents conclusions. 



2 A model of the growth of Saccharomyces cere- 
visiae 

Jones & Kompala [9] give the following model equations: 



'dt 



kLa[0 - O) - + X , 



dt ^ ^ ' V y2 Y. 

dei G /„ ^ \ ^ ^ (2.1) 



dt + G 

dea E 
au2- 



dt K2 + E 

deg G 



dt K3 + G 



dC 
'dt 



riVi + /3 j ei + a* 
riVi + 63 + a* , 



where X, G, E and O represent the concentrations of cell mass, glucose, ethanol (in 
gL~^) and dissolved oxygen (in mgL~^) in the culture volume, respectively. Each Cj 
represents the intracellular mass fraction of a key enzyme in the i*^ metabolic path- 
way. C represents the intracellular carbohydrate mass fraction. The subscripts 1,2 
and 3 correspond to the three pathways: fermentation, ethanol oxidation and glucose 
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oxidation, respectively; r, represents the growth rate on each pathway. Formulas for 
the growth rates and other functions are given in Appendix [Bl Details concerning 
the derivation of the model are found in [5] and expanded in the M.S. thesis of Jones 
[25]. 

In particular, the equation 

Vi = 7^ T , (2.2) 

max(ri, r2, rsj 

is introduced to model the regulation of enzyme activity by numerous biochemical 
mechanisms. Each Vi is a smooth function except at points where the two largest r, 
are equal; at these points each Vi has discontinuous derivatives with respect to some 
of the variables. All eight differential equations have at least one term containing a 
Vi and thus display the same lack of smoothness properties. Therefore the model is 
a piecewise-smooth, continuous system. Our goal is to understand the effects that 
this nonsmoothness has on the resulting dynamics. 

This paper focuses on variations in values of two parameters, namely the dilution 
rate, D (in h~^), and the dissolved oxygen mass transfer coefficient, fc^a (in h~^). 
Values for all other parameters are given in Appendix [B] and were kept constant 
throughout the investigations. 

A key property of the system (12.11) . is that the positive hyper-octant is forward 
invariant. That is, if the values of all variables are initially positive, they will remain 
positive for all time. This is, of course, a property that is required for any sensible 
model, since negative values of the variables are not physical. Furthermore, within 
the positive hyper-octant all trajectories are bounded forward in time. In other 
words, solutions always approach some attracting set which could be an equilibrium, 
a periodic orbit or a more complicated and possibly chaotic attractor. 



3 A bifurcation set 

In this section we describe a numerically computed bifurcation set for the system 
(12. ip . see Fig. [H We find a single, physically meaningful equilibrium (steady-state) 
except in small windows of parameter space between saddle-node bifurcations. For 
small values of D, this equilibrium is stable. If we fix the value of kia and increase 
D, the first bifurcation encountered is an Andronov-Hopf bifurcation (labelled HBi 
in Fig. [1]). Slightly to the right of HBi the equilibrium is unstable and solutions 
approach a periodic orbit or complicated attractor (see §1]). As the value of D 
is increased further a second Hopf bifurcation (HB2a or HB2b) is encountered that 
restores stability to the equilibrium (unless kia ^ 230, see below). 

Recall, the system is piecewise-smooth, continuous as a result of a maximum 
function in the rate coefficient, (12.21) . Since the switching manifold is codimension- 
one, it is a codimension-one phenomenon for the equilibrium to lie precisely on a 
switching manifold. Though an analytical formula for the equilibrium seems difficult 
to obtain, we have been able to numerically compute curves in two-dimensional 
parameter space along which this codimension-one situation occurs - the black curves 
in Fig. [H We refer to these as curves of discontinuity. The curves of discontinuity 
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Figure 1: A bifurcation set for the system (12. ip . HB - Hopf bifurcation, SN - saddle- 
node bifurcation. The parameter space is divided into three regions within each of 
which a different metabohc pathway is preferred at equihbrium. 



divide parameter space into three regions where one of the rj is larger than the other 
two, at the equilibrium. They may also correspond to discontinuous bifurcations, 
as described below. Physically, crossing a curve of discontinuity corresponds to a 
change in the preferred metabolic pathway at equilibrium. 

Fig. [2] shows an enlargement of Fig. [1] near two points on a curve of disconti- 
nuity. In panel A, along the curve of discontinuity, below the point (a) and above 
the point (c), there is no bifurcation. Between (a) and (c), numerically we have 
observed that a periodic orbit is created when the equilibrium crosses the switching 
manifold. Between (a) and (b), the orbit is unstable and emanates to the right of 
the curve of discontinuity. Between (b) and (c) the orbit is stable and emanates to 
the left. We refer to these bifurcations as subcritical and supercritical discontinuous 
Hopf bifurcations, respectively. The codimension-two point (b), is akin to a Hopf 
bifurcation at which the constant determining criticality vanishes. We expect that 
a locus of saddle- node bifurcations will emanate from this point, in manner similar 
to that in smooth systems. 

Two of the loci of smooth Hopf bifurcations, HB2a and HB2;,, collide with the 
curve of discontinuity at (a) and (c). Near these points these Hopf bifurcations 
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Figure 2: Magnified views of Fig. [H HB - Hopf bifurcation, SN - saddle-node bifur- 
cation. The curves of discontinuity are labelled by their corresponding bifurcations: 
NB - no bifurcation, DHB - discontinuous Hopf bifurcation (with criticality indi- 
cated), DSN - discontinuous saddle- node bifurcation. The dotted curves correspond 
to the grazing of a Hopf cycle with a switching manifold. There are three such 
curves, these emanate from the points (a), (c) and (d) (note, the curve that em- 
anates from (a) is barely distinguishable from HB2fe). In panel A, one point on a 
locus of saddle-node bifurcations of a Hopf cycle is shown with a red cross. 

are subcritical. Unstable periodic orbits emanate from the Hopf bifurcations and 
are initially of sufficiently small amplitude to not intersect a switching manifold. 
However, as we move away from the Hopf bifurcations, the amplitude of the Hopf 
cycles grow and they graze the switching manifold along the dotted curves in Fig. O 
No bifurcation occurs at the grazing because the system is continuous [26]. As we 
will show in Theorem 15. 4[ the grazing curves intersect the Hopf loci tangentially. 
The unstable cycles persist beyond grazing until they collide with a stable cycle in 
a saddle-node bifurcation. Loci of saddle-node bifurcations of periodic orbits are 
not shown in the figures because we have not been able to accurately numerically 
compute more than a single point (when kia = 240, see Fig. [2]A) on the curves 
due to the stiffness, non-smoothness and high dimensionality of the system fl2.ip . 
We expect one such curve to emanate from (c) and lie extremely close to the upper 
grazing curve as has recently been shown for two-dimensional systems |27j . 

Fig. [2]B, shows a second magnification of Fig. [T]near D = 0.225 and /c^a = 320. 
Loci of Hopf bifurcations and saddle-node bifurcations of the equilibrium have end- 
points at (d) and (e) that lie on a curve of discontinuity. We will show in ^ that 
bifurcations and dynamical behavior in neighborhoods of (d) and (e) are predicted 
by Theorems 15.41 and 15. respectively. To the left of the point (d), no bifurcation 
occurs along the curve of discontinuity. To the right of (e), points on the curve of 
discontinuity act as saddle-node bifurcations, hence we refer to these as discontin- 
uous saddle-node bifurcations. From the point (d) to very close to (e), the curve 
of discontinuity corresponds to a supercritical discontinuous Hopf bifurcation. A 
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Takens-Bogdanov bifurcation occurs at (f) where the Hopf locus, HB3, terminates 
at the saddle-node locus, and the point (g) corresponds to a cusp bifurcation [28] . 
Near the points (e) and (f) we believe there are a variety of additional bifurcations 
that we have not yet identified. The system exhibits stable oscillations in the re- 
gion between the smooth Hopf bifurcations and discontinuous Hopf bifurcations, 
but to our knowledge, oscillations in this parameter region have not been observed 
exp eriment ally. 



4 Simple and complicated stable oscillations 

The experimentally observed oscillations correspond to the region between HBi and 
HB2 in Fig. [H In this section, we will discuss the dynamics in this region in more 
detail. Fig. [3l shows a one parameter bifurcation diagram of the system (12.11) when 
k^a = 150. The black [magenta] dots represent local maxima [minima] O on the 
stable oscillating cycle. For values of D between about 0.0927 and 0.1172, all three 




Figure 3: A bifurcation diagram of the system (12. ip when kia = 150. The equilib- 
rium is colored blue when stable and red otherwise. Black [magenta] dots correspond 
to local maxima [minima] O that stable oscillations obtain. Two Hopf bifurcations 
are indicated by circles. The asterisk corresponds to where ri < r2 = at equilib- 
rium. The period of these oscillations is also indicated. 
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pathways are at some time preferred over one period of the solution, see Fig. |H In 
particular, very soon after the preferred pathway changes from glucose oxidation 
to fermentation (green to cyan in Fig. H]), the concentration of dissolved oxygen 
rebounds slightly before continuing to decrease. Thus local maxima appear below 
the equilibrium value in Fig. [31 For larger values of still to the left of the rightmost 
Hopf bifurcation, fermentation is no longer a preferred pathway at any point on the 
stable solution, and the lower local maximum is lost. Also, the absolute maximum 
undergoes two cusp catastrophes at ~ 0.111, 0.117. 
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Figure 4: Time series after transients have decayed when /c^a = 150, D = 0.1 in 
panel A and D = 0.12 in panel B. The solution is colored cyan when ri = r^aax, 
magenta when r2 = fmax and green when = rmax 



Different values of k^a yield similar bifurcation diagrams, we show a collection 
in Fig. [5l As a general rule there is a rapid change from a stable equilibrium to a 
large amplitude orbit near the leftmost Hopf bifurcation and as D is increased the 
amplitude and period of the orbit decrease. 

The behavior near the leftmost Hopf bifurcation is actually quite complex, as 
indicated in Fig. [6l which is a magnification of Fig. [31 Here the Hopf bifurcation 
is supercritical giving rise to a stable orbit which then undergoes a period-doubling 
cascade to chaos over an extremely small interval. The first period-doubling oc- 
curs at D ^ 0.092697 and the solution appears chaotic hj D = 0.092700. At 
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0.15 / 100 
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Figure 5: Bifurcation diagrams of (12.11) at six different values of k^a with the same 
color scheme as Fig. [3l A curve of discontinuity and Hopf loci shown in Fig. [1] are 
also included. 



D ^ 0.092701 the attractor suddenly explodes in size and the oscillation amplitude 
grows considerably. As D decreases toward this point we observe a period-adding 
sequence. Period- adding sequences are characterized by successive jumps in the 
period in a manner that forms an approximately arithmetic sequence. Such se- 
quences have been observed in models of many physical systems EHl EDI EI]- To 
our knowledge period-adding sequences are not completely understood, but seem to 
arise when periodic solutions interact with an invariant manifold of a saddle-type 
equilibrium giving rise to a Poincare map that is piecewise-smooth and often dis- 
continuous. Period-adding in one-dimensional piecewise-smooth maps has been the 
subject of recent research [201 E21 E3]- Dynamical behavior between period-adding 
windows (intervals of the bifurcation parameter within with the period undergoes no 
sudden change) is determined by the types and order of various local bifurcations. 
In Fig. [6], the period appears to go to infinity in the period-adding sequence. Within 
the extremely small regions between windows, we have identified period-doubling bi- 
furcations and complicated attractors although these attractors deviate only slightly 
from the observed periodic orbits. 
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Figure 6: A magnification of Fig. [3l near tlie leftmost Hopf bifurcation. Local 
minimums are not shown. 

5 Codimension-two discontinuous bifurcations 

This section studies dynamics near two of the codimension-two, discontinuous bi- 
furcation scenarios that were identified in ^ Adopting a general viewpoint, we will 
first unfold the simultaneous occurrence of a saddle-node bifurcation and a discon- 
tinuous bifurcation; our results are summarized in Fig. [7]\. The tangency illustrated 
in this figure matches our numerically computed bifurcation set, specifically point 
(e) of Fig. [2]B. Secondly we will unfold the simultaneous occurrence of a Hopf bifur- 
cation and a discontinuous bifurcation, see Fig. 03. This theoretical prediction also 
matches numerical results, specifically the points (a), (c) and (d) of Fig. [2l 

The results of this section are presented formally in Theorems 15.11 and 15 . 41 proofs 
of which are given in Appendix |Al Proofs to Lemmas 15.21 and 15.31 are described in 
[27] . Throughout this section we use arbitrary parameters and ri that do not relate 
to particular parameters of (12. ip . To clarify our notation, 0{i) \o{i)\ denotes terms 
that are order % or greater [greater than order z] in every variable and parameter on 
which the given expression depends. 

In the neighborhood of a single switching manifold an A^-dimensional, piecewise- 
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Figure 7: Unfoldings predicted by Theorems 15.11 and 15.41 SN - saddle-node bifur- 
cation, HB - Hopf bifurcation. Along the curve labelled "grazing" , the Hopf cycle 
intersects the switching manifold at one point. Along the T^-axis, an equilibrium lies 
on the switching manifold. 

smooth, continuous system may be written as 



where f^^\f^^^ : x ^ are and we assume that i/ : x M is 

sufficiently smooth (at least C^). The switching manifold is the parameter dependent 
set, S^^rj = {x E I H{x; /i, rj) = 0}. If, when (x; fi, rj) = (0; 0, 0), we have H = 
and VH ^ then locally iSo,o is a (A^ — l)-dimensional manifold intersecting the 
origin. Via coordinate transformations in a similar manner to those given in [26], 
we may assume, to order C^, that H is simply equal to x\. The higher order terms 
in H do not affect our analysis below; thus for simplicity, in what follows we will 
assume H is identically equal to x\. The switching manifold is then the plane xi = 
and we will refer to /'•■^^ as the left-half-system and f^^^ as the right-half-system. 

We assume that when /i = 77 = 0, the origin is an equilibrium. Since the origin 
lies on the switching manifold and (15.1 p is continuous, it is an equilibrium of both 
left and right-half-systems. We will assume 



that is, for the right-half-system, zero is not an associated eigenvalue of the origin 
when fi = rj = 0. Then by the implicit function theorem the right-half-system has an 
equilibrium, x*^^^ (//, ?]), with x*'''^^(0,0) = and that depends upon the parameters 
as a function in some neighborhood of the origin. As is generically the case, 
we may assume the distance the equilibrium is from the switching manifold varies 
linearly with some linear combination of the parameters. Without loss of generality 



X 



( /(^)(x;/i,?7), H{x;n,r]) 
1 /^^^a;;^,?]), H{x;n,r]) 



< 
> 



(5.1) 



det(D,/(^)(0;0,0))7^0; 



(5.2) 
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we may assume fi is a suitable choice; that is 



^-■'"''"■°>^0. (5,3) 

In this case, the imphcit function theorem imphes there is a function 0i : M ^ M 
such that x*^^\(f)i{r]) , r]) = 0. In other words when /x = the equilibrium lies 

on the switching manifold. By performing the nonlinear change of coordinates 

X ^ X — x*^^\(f)i{r]) , r]) , 
we may factor /i out of the constant term in the system, i.e. 

/(^)(0; ^, v) = f^^^Ko; v) = v) + oik) , 

where b is C^~^. Notice the transformation does not alter the switching manifold. 
The system is now 

f^\x;ii,ri), xi<0 



with 

f^\x; 12, 7]) = fibifi, ri) + Aifi, ri)x + Oi\x\^) + o(fc) , (5.5) 

where and Afi are N x N matrices that are C^~^ functions of /i and rj. 

Since (15. 4p is continuous, the matrices Ai and A^ have matching elements in all 
but possibly their first columns. It directly follows that the adjugate matrices (if A 
is non-singular, then adj(^) = det(74)74~^) of Al and Ar share the same first row 

^^ = e]"adj(^z.) = e]"adj(A^). (5.6) 

To understand the role of the vector, ^, consider equilibria of (15.41) . When ^ = 
1] = 0, the origin is an equilibrium. For small non-zero /x each half-system has an 
equilibrium, x*^^\ with first component 



det(y4j 



fi + 0{2) 

At=7?=0 



provided that det(Aj(0,0)) 7^ 0. Notice that our non-degeneracy assumption (15. 3p . 
is satisfied if ^"""(0,0)6(0,0) 7^ 0. If x*i^^ > 0, then x*^^^ is an equilibrium of the 
piecewise-smooth system (15.41) and is said to be admissible, otherwise it is virtual. 
Similarly x*^^^ is admissible if and only if x]^^^ < 0. Finally, notice if det{AL{0,0)) 
and det{Aji(0,0)) are of the same sign, then x*^^^ and x*^^^ are admissible for dif- 
ferent signs of fi, whereas if det(^L(0,0)) and det{Aji{0,0)) have opposite sign, 
x*(^) and x*^^^ are admissible for the same sign of n. The former case is known as 
persistence; the latter is known as a non-smooth fold |20j . 

The following theorem describes dynamical phenomena near = r] = when 
det(AL(0,0)) = 0. 
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Theorem 5.1 

Consider the system ( [5.^[ j with ( (5.51) and assume that N > 1 and k > A. Suppose 
that near {fi,r]) = (0,0), A^dj,,?]) has an eigenvalue X{fj,,ri) G M with the associated 
eigenvector, v{^,ri). In addition, suppose 

i) A(0,0) = zs of algebraic multiplicity 1 and is the only eigenvalue 0/^4^(0,0) 
with zero real part. 

Then f (0, 0) has a non-zero first component and the magnitude of v(0, 0) may be 
scaled such that, ^(0, 0)"^v(0, 0) = 1. Finally suppose that 



1(0,0,0) 



ii) f (0,0)^0, 

m; a, = e{{Dlf^''^){v.,v))\ 

iv) ^(0,0)^6(0,0) ^ 0. 
Then there exists a unique C^'"^ function h : 



with h{0) = h'{0) = and 



^ dq ' 



2aoCb 



(5.7) 



(0,0) 



such that in a neighborhood of {fi,r]) = (0,0), the curve /i = /i(?7) corresponds to a 
locus of saddle-node bifurcations of equilibria of ( [5.^[ j that are admissible when 



sgn(?7) = sgn(aot;iff 



dn 1 1 (0,0) 



(5i 



A proof of Theorem 15.11 is given in Appendix |Al The theorem implies a bifurca- 
tion diagram hke that depicted in Fig. [7|A. In particular the curve of saddle-node 
bifurcations is tangent to the T^-axis as shown. 

The second theorem. Theorem 15.41 describes dynamical phenomena near = 
1] = when ^4^(0,0) has a purely imaginary complex eigenvalue pair. The method 
of proof is essentially a standard dimension reduction by restriction to the center 
manifold of the left-half-system The dynamics of the resulting planar system are 
determined by the following two lemmas. Lemma 15.31 provides a transformation 
of the planar system to observer canonical form [20j. Lemma 15.21 describes local 
dynamics of the planar system in this canonical form. 



Lemma 5.2 

Consider the two-dimensional (k > 5) system 



X 









f{x,y;n,r)) 



Suppose, 

i) (5(0,0) =uj^ foru > 0, 








+ 









r] 1 
-6{fi,v) 





X 







+ 0{\x,y\') + o{k) 



(5.9) 
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ii) ao 7^ 0, where 

f^O ^ fxxx ~l~ Qxxy ~\~ ^ fxyy~\~^ 9yyy ) •^^y ( ~^ -A/y ) 

+ '77:9xy{ nSxx + + K-f^^Sxx ~ '-^ fyydyy) (5.10) 

evaluated at {x,y; fj,,?]) = (0,0; 0,0). Then, near {x,y; n,!]) = (0,0; 0,0), there is a 
unique equilibrium given by functions 

x*{fi,ri) = --^^ + 0(2), 
y*ifi,v) = 0(2), 
and there exist C^~^ , C^~^ functions /ii, /12 : M ^ M respectively, with 



hlifJ') = -^{fxx+9xy) 



UJ 

2 I /n/',,3 



At + 0(^'), (5.11) 

(0,0;0,0) 



hM = /ii(/i)-^/i^ + 0(/i^), (5.12) 



such that for small fi, the curve rj = hi{^) corresponds to Andronov-Hopf bifurcations 
of (x*, y*y and the curve rj = /i2(/^) corresponds to associated Hopf cycles intersect- 
ing the y-axis tangentially at one point. The Hopf bifurcations are supercritical if 
Co < and subcritical if > 0. The Hopf cycle lies entirely in the left-half-plane if 
and only if fi > and 77 lies between hi{fi) and ^2(/^)- 

Lemma 5.3 

Consider the system ( [5.^| ) with Ii5.5\) and assume that N = 2 and k > 5. Suppose 
Al{^, rf) has complex eigenvalues A± = z/ ± io; with 

i) zy(0,0) = anc? cu(0,0) > 0, 

ii) |(0,0)^0, 

ill) ^^(0, 0)6(0,0) ^ 0. 

Then there is a nonlinear transformation not altering the switching manifold such 
that the left- half -flow of ( [5./^[ j is given by Ii5. 9\) . All conditions in Lemma lSTB will be 
satisfied except possibly the non- degeneracy condition, oq 7^ 0. 

Theorem 5.4 

Consider the system ( [5.^[ ) with Ii5.5\) and assume that N > 2 and k > 6. Suppose 
near {^,ri) = (0,0), AL{fi,r]) has complex eigenvalues X± = u±iuj with associated 
eigenvectors, z± = u'^^^ ± iu^'^\ Suppose 

i) z^(0,0) = 0, u;(0,0) > 0, and Al{0,0) has no other eigenvalues on the imagi- 
nary axis, 

11) |(0,0)^0, 
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m) ^^(0, 0)6(0,0) ^ 0, 



iv) either ul or ul is non-zero. 

Then, in the extended coordinate system {x, fi, rj), there exists a C^~^ four-dimensional 
center manifold, W^, for the left-half- system that passes through the origin and is 
not tangent to the switching manifold at this point. Furthermore, 3i 1, such that 
in the coordinate system 



Xi 




Xi 


X2 




Xi 



the left-half-fiow of ( [5.^[ ) restricted to is given by / (5.g|) (with "hatted" variables) 
and all conditions in Lemma 15.51 will be satisfied. 

See Appendix |A] for a proof. Theorem 15.41 implies that a n-dimensional system near 
a codimension-two point with a simultaneous Hopf and discontinuous bifurcation 
will have a bifurcation diagram like that shown in Fig. [7|B- In particular, the curves 
of grazing and Hopf bifurcations are tangent to one another at the origin. For this 
scenario in two-dimensions it is known that a curve of saddle-node bifurcations of 
the Hopf cycle may exist very close to the grazing curve [27] . We have not been able 
to extend this result to higher dimensions. 

6 Conclusions 

In this paper we investigated the onset of stable oscillations and more complex 
behavior in a model of S. cerevisiae growth taken from [9J. The model assumes 
an instantaneous switching between competing metabolic pathways resulting in a 
piecewise-smooth, continuous system of ODEs. In this paper we identified a variety 
of discontinuity induced bifurcations. 

The model exhibits stable oscillations that arise from Andronov-Hopf bifurca- 
tions for intermediate values of the dilution rate, D; these have also been observed 
experimentally. As D grows, the oscillation amplitude suddenly jumps to a much 
larger value just slightly beyond the Hopf bifurcation. We do not have a detailed 
explanation for this sudden amplitude change. As D is increased further, the result- 
ing stable orbits undergo a complex sequence of bifurcations causing their periods 
and amplitudes to decrease, until a second Hopf bifurcation occurs resulting again 
in a stable equilibrium. 

For the model (12.11) . a change in the preferred metabolic pathway at an equi- 
librium results in the loss of differentiability for orbits in its neighborhood. The 
result is often a discontinuity induced bifurcation, and we have identified discontin- 
uous saddle-node and Hopf bifurcations. The system also exhibits codimension-two 
bifurcations that correspond to simultaneous discontinuous saddle-node and Hopf 
bifurcations. We have provided a rigorous unfolding of these scenarios from a general 
viewpoint. 

While the behavior that we studied are specific to piecewise smooth, continuous 
models, a model in which this simplification is relaxed should still have much the 
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same behavior. For example if the relative strength of two competing pathways 
reverses, then exponential growth will lead to the dominance of one over the other 
over a small parameter range and a short timescale. Though the bifurcations generic 
to smooth systems are restricted relative to those of discontinuous systems, a rapid 
sequence of these bifurcations over a small range of parameters may lead to the same 
behavior on a rougher scale as the discontinuous one. This can be seen, for example, 
in the simplest models such as a smoothed one- dimensional tent map. 



A Proofs 



Proof of Theorem 15.11 

Recall, for any N x N matrix. A, 

adj(A)A = det(A)J . 
By putting A = Al{0, 0), multiplying on the left by ej and using (15. 6p we obtain 

e{0,0)AL{0,0) = 0. 

Thus ^"""(0, 0) is the left eigenvector of Al{0, 0) for A(0, 0) = 0. Consequently we may 
indeed choose the length of the right eigenvector f (0, 0) such that ^"""(0, 0)^;(0, 0) = 1. 

Now we show v{0, 0) has a non-zero first element. Suppose for a contradiction, 
fi(0, 0) = 0. Let Bij denote the (A^— 1) x [N — 1) matrix formed by removing the i^^ 
row and column from Al(0,0). Let v = {v2{0,0) . . . Vn{0,Q)V G ffi^"^ Then 
V ^ and Bnv = for each i. Therefore each det{Bii) = 0, i.e., each element in 
the first column of the cof actor matrix of ^4^(0, 0) is zero. Thus ,^"^(0, 0) = which 
is a contradiction. Therefore 

^i(0,0)^0. 

Let v^^\ . . . , v^^^ be N generalized eigenvectors of Al(0, 0) that form a basis of 



with = v(0,0). 
coordinates 

Let 



Let V 



X 



] . We introduce the linear change of 
X . (A.l) 



X 













V-'f^iVx;fi,v) 





(A.2) 



denote the (A^ + 2)-dimensional C'^ extended left-half- flow in the basis of generalized 
eigenvectors. The Jacobian 



DF(0;0,0) 



y-Mi(o,o)y 

















(A.3) 
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has a three-dimensional nuUspace, E'^. The {N + 2)-dimensional vectors 



{ni,n2,ns} = {ei, 



1 




where 77-2 is a generahzed eigenvector and G 
manifold, W^, is tangent to i?^, thus on W^, 



) eN+2} 1 

M^, span E"^. The local center 



X = H{xi] /i, rf) = xiCi - fiC + 0(2) , 

where C G is equal to Lp except that its first element is zero. Notice ^"""(0, 0)V = ej 
thus 

= ejx = ^{0, 0)Vx = ^{0, 0)x . 
Restricted to the dynamics (\A.2\i become the C'^^^ system 



£1 = tie{0,0)b{fi,v) + e{0,0)ALifi,T^)VH{x,;^,r]) 
+ e{0,0)9''{VH{xufi,rj);fi,r,) , 



(A.4) 



where g^{x; fi, rj) denotes all terms of that are nonlinear in x. By expanding each 
term in (1A.4P to second order we obtain 



xi{xi] jj, r]) = ci/x + C2xl + csXiT] + C4X1H + + cefiT] + 0(3) 



(A.5) 



where, in particular 



ci = e(0,0)"^6(0,0) 
«o 
2 ' 
dX 



C3 



dr] 



(0,0). 



Let xl be an equihbrium of (lA.Sp . Since Ci 7^ by hypothesis, by the implicit 
function theorem there exists a unique C''~^ function, /Ueq(xJ,?]) such that Xi(f^; 
^eqixl,!]),!]) = and 

fleaix*!, V) = -—xf - —xIt] + 0(3) . 
Ci Ci 

Near {n,r)) = (0,0), the linearization about the equilibrium xl, has an associated 
eigenvalue of exactly when 







dxi 
dxi 



[X 



t; fJ-eqixl, V),V) = 2C2Xl + CsT] + 0(2) . 



(A.6) 



Since ao 7^ by hypothesis, the implicit function theorem again implies there exists 
a unique C^~^ function, /i : R — >• R such that, (IA.6P is satisfied when xl = h{rj) and 



Hv) 



_C3_ 
'2C2 
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(A.7) 



The function 



4ciC2 

is C^~'^. We now show saddle-node bifurcations occur for the left-half- flow on the 
curve = h{rj) when r] is small by verifying the three conditions of the saddle-node 
bifurcation theorem, see for instance i34i: 



i) by construction, Dxf^{x*;h{r]),r]) has a zero eigenvalue of algebraic multi- 
plicity 1, and there are no other eigenvalues with zero real part when r] is 
sufficiently small, 

ii) w{fi,r])^{x*;h{7]),ri) = ^^(o, 0)6(0, 0) + 0{r]) ^ where w{^,ri) is the left 
eigenvector of AL{fi,r]) for A(/i,r^), 

iii) w{fi, r]){Dlf^{x*] h{r]), r]){v{fi, v))) = ao + 0{r]) 0. 
Finally, notice that on 

xi = eJVH{xu /i, rj) = vi{0, 0)xi + 0(2) ^ , 
and by (]A.7p . when /i = h{r]) 

C3Ui(0,0) 
2C2 



Thus the equilibrium at the saddle- node bifurcation is admissible when sgn(77) 

dr) ' 



sgn(c2C3z;i(0,0)) = sgn(aofi^^ 



(0,0) 

□ 

Proof of Theorem 15.41 

Since the real and imaginary parts of the eigenvectors u^^^ and -u*-^-*, are linearly 
independent, there exists an z 7^ 1 such that 



U 



(1) (2) 

u\ u\ 

(1) (2) 



is non-singular. For the remainder of this proof we will set z = 2, w.l.o.g. Define 
two new vectors by 

[^(1) = [^(1) ^(2)] ^ (A.9) 

let 

and introduce the new coordinate system 

X = V-^x . (A.IO) 

Note the inclusion of the matrix, U"^, in f lA.Op allows for simplification below, in 
particular, eJV = eJV^^ = ej for i = 1,2. The (A^ + 2)-dimensional C'' extended 
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left-half-flow in the new coordinates is given by (1A.2I) as before. The Jacobian, 



DF(0; 0,0), (1A.3I1 . has a four- dimensional linear center manifold, E^, spanned by 



{ni,n2,n3,ni} = {ei,e2. 



-{AL{0,0)V)-'b 
1 




eN+2} ■ 



Notice is not tangent to the switching manifold by condition (iv) of the theorem. 
On the local center manifold 

X = H{xi, X2] /i, T]) = xiCi + £262 - Ca* + 0{2) , 

where C ^ is equal to (^^(0, 0)V)~^b except that its first two elements are zero. 
The dynamics on are described by 



Xi 




\ej] 


( 


. ^2 . 




.el _ 





V-^^b+V~^ALVH{xi,X2]^i,r]) 
+ \/-V^)(l^//(xi,X2;/i,?7);/i,?7) ) , 



where g^^'' represents all terms of /'•^^ that are nonlinear in x. By using (1A.9P and 



where 

we obtain 

where 

and 
h = 



D 



Xx 
X2 

At. 



Xi 

£2 



0(2) 



(A.ll) 



el 



V-^Ajy\ex 62] = UDU-^ , 



eJ 
el 



Al 0---0 

Finally, it is easily verified that 



V~'AlVC 



V-'Al [v'^'^ ^;(2)0---0] Al'b 



Atb. 



where ^ = eJad}{AL). 
□ 
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B Functions and Parameter Values 



The following is a complete list of functions that are present in the model 

fJ'i — A*i,max . ^ j ^ ^j'^t'^ i 

a + a* 
G 



E O 

GO 

- ^'''Ks + GKo, + 0' 

^ = 1,2,3, 
z= 1,2,3 



r,; 



maxj Tj 



Go 


10 gL-i 


Ki 


0.05 gL-i 




0.16 gg-i 


K2 


0.01 gL-i 




0.75 gg-i 


Ks 


0.001 gL-^ 




0.60 gg-^ 


Ko, 


0.01 mgL"^ 


01 


0.403 


Ko, 


2.2 mgL"^ 


02 


2000 mgg-i 


71 


10 


03 


1000 mgg-i 


72 


10 


04 


0.95 


73 


0.8 gg-i 


0* 


7.5 mgL~^ 


A''l,max 


0.44 h-i 


a 


0.3 gg-^h-^ 


/-^2,max 


0.19 h-i 


a* 


0.1 gg-^h-^ 


A''3,max 


0.36 h-i 


P 


0.7 h-i 







Table 1: Parameter values used for all numerical investigations in this paper. 
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